## Created by EmpowerStats @ Mon, 17 Feb 25 21:05:45 +0800## 
#******************** Regarding ALL Following R Code   ******************** 
#*** COPYRIGHT (c) 2010, 2021 X&Y Solutions, ALL RIGHT RESERVED *********** 
#******************* www.EmpowerStats.com ********************************* 
#************************************************************************** 
Sys.setlocale(category = 'LC_ALL', locale = 'English_United States.1252'); 
.libPaths(file.path(R.home(),'library')); 
library(doBy); 
options(timeout=600); 
library(plotrix); 
library(stringi); 
library(stringr); 
library(survival); 
library(rms); 
library(nnet); 
library(car); 
library(mgcv); 
pdfwd<-6; pdfht<-6; 
load('/Users/liukangning/Desktop/EmpowerROS/Analysis/cd/datamergeall.Rdata'); 
if (length(which(ls()=='EmpowerStatsR'))==0) EmpowerStatsR<-get(ls()[1]); 
names(EmpowerStatsR)<-toupper(names(EmpowerStatsR)); 
originalVNAME<-names(EmpowerStatsR); 
ofname<-'cd_10_tbl'; 
vname<-c(NA,'EDUCATION','EDUCATION.1','EDUCATION.2','EDUCATION.3','EDUCATION.4','EDUCATION.5','DPQ','AFT','CERAD','DSST')[-1]; 
vlabel<-c(NA,'EDUCATION','  1','  2','  3','  4','  5','DPQ','AFT','CERAD','DSST')[-1]; 
varused4this <- c('EDUCATION','DPQ','AFT','CERAD','DSST'); 
pkgs<-c('mgcv','geepack','gdata'); 
for (g in pkgs) {  
if (!(g %in% rownames(installed.packages()))) install.packages(g,repos='https://cloud.r-project.org'); 
}
library(mgcv); 
library(geepack); 
library(gdata); 
WD <- EmpowerStatsR; rm(EmpowerStatsR); gc(); 
title<-'单因素分析'; 
wd.subset=''; 
svy<- 0; 
weights.var<- NA; 
yvname<-c('DSST','CERAD','AFT'); 
ydist<-c('gaussian','gaussian','gaussian'); 
ylink<-c('identity','identity','identity'); 
ylv<-c(0,0,0); 
xvname<-c('DPQ'); 
sxf<-c(0); 
xlv<-c(0); 
subjvname<- NA; 
gee.TYPE<-NA; 
cox<- 0; 
timevar<- NA; 
vname.start<- NA; 
avname<-c(); 
bvar<-'EDUCATION'; blv<- 5; 
colvname<- NA; 
prn<-1; 
dec<-1; 
##R package## mgcv geepack gdata ##R package##;
pvformat<-function(p,dec) {
  pp <- sprintf(paste("%.",dec,"f",sep=""),as.numeric(p))
  if (is.matrix(p)) {pp<-matrix(pp, nrow=nrow(p)); colnames(pp)<-colnames(p);rownames(pp)<-rownames(p);}
  lw <- paste("<",substr("0.00000000000",1,dec+1),"1",sep="");
  pp[as.numeric(p)<(1/10^dec)]<-lw
  return(pp)
}
numfmt<-function(p,dec) {
  if (is.list(p)) p<-as.matrix(p)
  pp <- sprintf(paste("%.",dec,"f",sep=""),as.numeric(p))
  if (is.matrix(p)) {pp<-matrix(pp, nrow=nrow(p));colnames(pp)<-colnames(p);rownames(pp)<-rownames(p);}
  pp[as.numeric(p)>10000000]<- "inf."
  pp[is.na(p) | gsub(" ","",p)==""]<- ""
  pp[p=="-Inf"]<-"-Inf"
  pp[p=="Inf"]<-"Inf"
  return(pp)
}
varstats<-function(var,vlvl,dec) {
  if (length(vlvl)==1 & vlvl[1]==0) {
    return(paste(numfmt(mean(var,na.rm=TRUE),dec),numfmt(sd(var,na.rm=TRUE),dec),sep=" &#177 "))
  } else {
    a<-table(var)
    b<-matrix(paste(a, " (", numfmt(a/sum(a)*100,dec), "%)",sep=""),ncol=1)
    return(c(" ",b[match(vlvl,names(a))]))
  }
}
mat2htmltable<-function(mat) {
  t1<- apply(mat,1,function(z) paste(z,collapse="</td><td>"))
  t2<- paste("<tr><td>",t1,"</td></tr>")
  return(paste(t2,collapse=" "))
}
setgam<-function(fml,yi) {
  if (ydist[yi]=="") ydist[yi]<-"gaussian"
  if (ydist[yi]=="exact") ydist[yi]<-"binomial"
  if (ydist[yi]=="breslow") ydist[yi]<-"binomial"
  if (ydist[yi]=="gaussian") mdl<-try(gam(formula(fml),weights=wdtmp$weights,data=wdtmp, family=gaussian(link="identity")))
  if (ydist[yi]=="binomial") mdl<-try(gam(formula(fml),weights=wdtmp$weights,data=wdtmp, family=binomial(link="logit")))
  if (ydist[yi]=="poisson") mdl<-try(gam(formula(fml),weights=wdtmp$weights,data=wdtmp, family=poisson(link="log")))
  if (ydist[yi]=="gamma") mdl<-try(gam(formula(fml),weights=wdtmp$weights,data=wdtmp, family=Gamma(link="inverse")))
  if (ydist[yi]=="negbin") mdl<-try(gam(formula(fml),weights=wdtmp$weights,data=wdtmp, family=negbin(c(1,10), link="log")))
  return(mdl)
}
setgee<-function(fml,yi) {
  if (ydist[yi]=="") ydist[yi]<-"gaussian"
  if (ydist[yi]=="exact") ydist[yi]<-"binomial"
  if (ydist[yi]=="breslow") ydist[yi]<-"binomial"
  if (ydist[yi]=="gaussian") md<-try(geeglm(formula(fml),id=wdtmp[,subjvname],corstr=gee.TYPE,family="gaussian",weights=wdtmp$weights,data=wdtmp))
  if (ydist[yi]=="binomial") md<-try(geeglm(formula(fml),id=wdtmp[,subjvname],corstr=gee.TYPE,family="binomial",weights=wdtmp$weights,data=wdtmp))
  if (ydist[yi]=="poisson") md<-try(geeglm(formula(fml),id=wdtmp[,subjvname],corstr=gee.TYPE,family="poisson",weights=wdtmp$weights,data=wdtmp))
  if (ydist[yi]=="gamma") md<-try(geeglm(formula(fml),id=wdtmp[,subjvname],corstr=gee.TYPE,family="Gamma",weights=wdtmp$weights,data=wdtmp))
  if (ydist[yi]=="negbin") md<-try(geeglm.nb(formula(fml),id=wdtmp[,subjvname],corstr=gee.TYPE,weights=wdtmp$weights,data=wdtmp))
  return(md)
}
setglm<-function(fml,yi) {
  if (ydist[yi]=="") ydist[yi]<-"gaussian"
  if (ydist[yi]=="exact") ydist[yi]<-"binomial"
  if (ydist[yi]=="breslow") ydist[yi]<-"binomial"
  if (ydist[yi]=="gaussian") md<-try(glm(formula(fml),family="gaussian",weights=wdtmp$weights,data=wdtmp))
  if (ydist[yi]=="binomial") md<-try(glm(formula(fml),family="binomial",weights=wdtmp$weights,data=wdtmp))
  if (ydist[yi]=="poisson") md<-try(glm(formula(fml),family="poisson",weights=wdtmp$weights,data=wdtmp))
  if (ydist[yi]=="gamma") md<-try(glm(formula(fml),family="Gamma",weights=wdtmp$weights,data=wdtmp))
  if (ydist[yi]=="negbin") md<-try(glm.nb(formula(fml),weights=wdtmp$weights,data=wdtmp))
  return(md)
}
mdl2oo<-function(mdl, xxname, opt) {
  if (is.na(mdl[[1]][1])) return(list(rep("",times=length(xxname)),""))
  if (substr(mdl[[1]][1],1,5)=="Error") return(list(rep("",times=length(xxname)),""))
  gs<-summary(mdl); print(mdl$formula); print(gs)
  if (opt=="gam") {gsparm <- gs$p.table;tmpn<-gs$n;
  } else {gsparm <- gs$coefficients;tmpn <- sum(gs$df[c(1,2)]);}
  gsp<-gsparm[match(xxname,rownames(gsparm)),]
  if (length(xxname)==1) {beta<-gsp[1]; se<-gsp[2]; pv<-gsp[4]; 
  } else {beta<-gsp[,1]; se<-gsp[,2]; pv<-gsp[,4]; }
  ci1<- beta-1.96*se; ci2<- beta+1.96*se
  pvx<-substr(rep("****",length(pv)),1,(pv<=0.05)+(pv<=0.01)+(pv<=0.001)) 
  if (colprn==3) {pvv<-pvx;} else {pvv<-pvformat(pv,dec+2);}
  if ((colprn!=2) & (gs$family[[2]]=="log" | gs$family[[2]]=="logit")) {
    o1<-paste(numfmt(exp(beta),dec)," (",numfmt(exp(ci1),dec),", ",numfmt(exp(ci2),dec),")",sep="")
  } else {
    if (colprn<3) {o1<-paste(numfmt(beta,dec), " (",numfmt(ci1,dec),", ",numfmt(ci2,dec),")",sep="")
    } else {o1<-paste(numfmt(beta,dec), "+",numfmt(se,dec),sep="");}
  }
  o1<-paste(o1,pvv); o1[is.na(beta)]<-NA;
  if (length(xxname)>1) {
    o1[1]<-""; o1[2]<-"0";
    if (gs$family[[2]]=="log" | gs$family[[2]]=="logit") o1[2]<-"1.0"
  }
  return(list(o1,tmpn))
}
recodevar <- function (var,oldcode,newcode) {
  tmp.v <- var
  nc.tmp <- length(oldcode)
  for (i in (1:nc.tmp)) {tmp.v[(var==oldcode[i])]=newcode[i]}
  if (is.factor(tmp.v)) {tmp.v1<-as.numeric(as.character(tmp.v))} else {tmp.v1<-as.numeric(tmp.v)}
  rm(tmp.v);  return(tmp.v1)
} 
rankvar <- function(var, num) {
  qprobs <- 1/num
  if (num>2) {for (i in (2:(num-1))) {qprobs <- c(qprobs, 1/num * i) } } 
  outvar <- rep(0, times=length(var))
  outvar[is.na(var)] <- NA
  cutpoints <- quantile(var,probs=qprobs, na.rm=TRUE)
  for (k in (1:length(cutpoints))) { outvar[var>=cutpoints[k]] <- k; }
  return(outvar)
}
removeNA<-function(i,j,wdf) {
 vvv<-c(yvname[i],xvname[j],avname,subjvname,colvname,bvar,vname.start,timevar); 
 vvv<-vvv[!is.na(vvv)]; vvv<-vvv[vvv>" "]
 tmp<-is.na(wdf[,vvv]); 
 return(wdf[apply(tmp,1,sum)==0,])
}
if (!is.na(weights.var)) {weights<-WD[,weights.var];} else {weights<-1;}
WD<-cbind(WD,weights);
vlabelN<-(substr(vlabel,1,1)==" ");
vlabelZ<-vlabel[vlabelN];vlabelV<-vlabel[!vlabelN]
vnameV<-vname[!vlabelN];vnameZ<-vname[vlabelN]
w<-c("<!DOCTYPE html><html lang='zh'><head><meta charset='utf-8'></head><body>")
w<-c(w,paste("<h2>", title, "</h2>"))
if (length(avname)>0) { 
  if (sum((saf=="s" | saf=="S") & alv>0)>0) w<-c(w,"</br>Spline smoothing only applies for continuous variables")
  if (!is.na(subjvname) & (sum((saf=="s" | saf=="S") & alv==0)>0)) w<-c(w,"</br>Generalized estimate equation could not be used with spline smoothing terms")
}
allvname<-c(yvname,xvname,colvname,bvar,avname,subjvname,vname.start,timevar,"weights"); 
allvname<-allvname[!is.na(allvname)];
WD<-WD[,allvname];
if (!is.na(subjvname)) {
  if (length(avname)>0) saf<-rep(0,length(saf)); 
  WD<-WD[order(WD[,subjvname]),];
}
sxf<-as.numeric(sxf);sxf[is.na(sxf)]<-0;
if (sum(sxf>1 & xlv>0)>0) w<-c(w,"Categorizing only applies to continuous variables");
if (sum(sxf>1 & xlv==0)>0) {
     t.xname<-NA;t.xlv<-NA; nx<-length(xvname)
     for (i in 1:nx) {
       if (sxf[i]>1 & xlv[i]==0) {
          tmp.Xi<- rankvar(WD[,xvname[i]],sxf[i])
          tmp.newcode <- tapply(WD[,xvname[i]],tmp.Xi,function(z) median(z,na.rm=TRUE))
          tmp.low <- tapply(WD[,xvname[i]],tmp.Xi,function(z) min(z,na.rm=TRUE))
          tmp.upp <- tapply(WD[,xvname[i]],tmp.Xi,function(z) max(z,na.rm=TRUE))
          tmp.Xi2<- recodevar(tmp.Xi,(1:sxf[i])-1,tmp.newcode)
          tmp.Xi<-cbind(tmp.Xi,tmp.Xi2)
          tmp.NM<-paste(xvname[i],c("grp","grp.cont"),sep=".")
          colnames(tmp.Xi)<-tmp.NM
          WD<-cbind(WD,tmp.Xi)
          t.xname<-c(t.xname,tmp.NM)
          t.xlv<-c(t.xlv,sxf[i],0)
          vnameV<-c(vnameV,tmp.NM)
          vlabelV<-c(vlabelV,paste(vlabelV[vnameV==xvname[i]],c("group","group trend")))
          vnameZ<-c(vnameZ,paste(tmp.NM[1],(1:sxf[i])-1,sep="."))
          vlabelZ<-c(vlabelZ,paste(tmp.low,"-",tmp.upp))
       } else {
          t.xname<-c(t.xname,xvname[i]); t.xlv<-c(t.xlv,xlv[i])
       }
     }
     xvname<-t.xname[-1]; xlv<-t.xlv[-1];
}
fml0<-""; na=0; avb=""; smoothav<-0; 
if (length(avname)>0) {
  na<-length(avname)
  avb<-vlabelV[match(avname,vnameV)]; avb[is.na(avb)]<-avname[is.na(avb)]
  avname_ <- avname 
  smoothavi<-((saf=="s" | saf=="S") & alv==0)
  smoothav<-sum(smoothavi)
  smoothavname<-avname[smoothavi]
  avname_[smoothavi]<-paste("s(",avname[smoothavi],")",sep="")
  avb1<-avb
  avb1[smoothavi]<-paste(avb[smoothavi],"(Smooth)",sep="")
  avname_[alv>0]<-paste("factor(",avname[alv>0],")",sep="")
  fml0<-paste("+",paste(avname_,collapse="+"))
}
ny=length(yvname); nx=length(xvname); 
xb<-vlabelV[match(xvname,vnameV)]; xb[is.na(xb)]<-xvname[is.na(xb)]
yb<-vlabelV[match(yvname,vnameV)]; yb[is.na(yb)]<-yvname[is.na(yb)]
xvname_ <- xvname
xvname_[xlv>0]<-paste("factor(",xvname[xlv>0],")",sep="")
xxname_<-list(NA); xxlbl_<-list(NA); xxlvl_<-list(NA)
for (j in (1:nx)) {
  if (xlv[j]==0) {
    xxname_[[j+1]]<-xvname[j];xxlbl_[[j+1]]<-xb[j];xxlvl_[[j+1]]<-0
  } else {
    xxlvl_[[j+1]]<-levels(factor(WD[,xvname[j]]))
    tmp<-paste(xvname[j],".",xxlvl_[[j+1]],sep="")
    xxlbl_[[j+1]]<-c(xb[j],vlabelZ[match(tmp,vnameZ)])
    xxlbl_[[j+1]]<-paste(c("",rep("&nbsp&nbsp",length(xxlbl_[[j+1]])-1)),xxlbl_[[j+1]])
    xxname_[[j+1]]<-c(xvname[j],paste("factor(",xvname[j],")",xxlvl_[[j+1]],sep=""))
  }
}
xxname_<-xxname_[-1]; xxlbl_<-xxlbl_[-1]; xxlvl_<-xxlvl_[-1];
fml0<-paste(xvname_,fml0)
if (is.na(colvname)) {
  nclv<-1; clvb<-"Total"; clvb_<-"Total"
} else {
  clv<-levels(factor(WD[,colvname])); nclv<-length(clv)+1
  clvb_<-vlabelZ[match(paste(colvname,".",clv,sep=""),vnameZ)]; clvb_[is.na(clvb_)]<-clv[is.na(clvb_)];
  clvb<-c(paste(vlabelV[vnameV==colvname],clvb_,sep="="),"Total");
  clvb_<-c(clvb_,"Total")
  WD<-WD[!is.na(WD[,colvname]),]
} 

if (is.na(bvar)) {
  blvb<-"Total"; blvb_<-"Total"
} else {
  blv<-levels(factor(WD[,bvar])); nblv<-length(blv)+1
  blvb_<-vlabelZ[match(paste(bvar,".",blv,sep=""),vnameZ)]; blvb_[is.na(blvb_)]<-blv[is.na(blvb_)];
  blvb<-c(paste(vlabelV[vnameV==bvar],blvb_,sep="="),"Total");
  blvb_<-c(blvb_,"Total")
  WD<-WD[!is.na(WD[,bvar]),]
}
opt<-ifelse(!is.na(subjvname), "gee", ifelse(smoothav>0, "gam", "glm")) ;
colprn<-prn; if (is.na(colprn)) colprn<-1
outprn<-0; 
if (!is.na(bvar)) outprn<-1
if (ny>4 & ny>max(xlv) & nx==sum(xlv>1) & nx<3 & is.na(bvar)) outprn<-2
if (ny>4 & ny>nx & nx==sum(xlv==0) & nx<5 & is.na(bvar)) outprn<-3
sink(paste(ofname,".lst",sep=""))
if (outprn==0) {
  tt<-c(" ","Statistics",yb); nn<-c("","Exposure",yb); 
  for (k in (1:nclv)) {
    wdtmp0<-WD;
    if (!is.na(colvname)) {
      tt<-rbind(tt,c(clvb[k],rep(" ",ny+1)))
      if (k<nclv) wdtmp0<-WD[WD[,colvname]==clv[k],];
      print(paste("Stratified by",colvname, ":", clvb[k]))
    } 
    for (j in (1:nx)) {
        colj<-cbind(xxlbl_[[j]],varstats(wdtmp0[,xvname[j]],xxlvl_[[j]],dec))
        nnj <- c(clvb_[k],xb[j])
        for (i in (1:ny)) {
          wdtmp<-removeNA(i,j,wdtmp0)
          fml<-paste(yvname[i],"~",fml0[j]);
          if (k==nclv & !is.na(colvname)) fml<-paste(fml,"+factor(",colvname,")",sep="")
          if (opt=="gam") tmp.mdl<-setgam(fml,i)
          if (opt=="gee") tmp.mdl<-setgee(fml,i)
          if (opt=="glm") tmp.mdl<-setglm(fml,i)
          tmpooi<-mdl2oo(tmp.mdl,xxname_[[j]],opt)
          colj<-cbind(colj,tmpooi[[1]])
          nnj<-c(nnj,tmpooi[[2]])
        }
        tt<-rbind(tt,colj); nn<-rbind(nn,nnj)
    }
  }
  if (is.na(colvname)) {nn<-nn[,-1];} else {nn[1,1]<-vlabelV[vnameV==colvname];}
  w<-c(w,"<table border=3>", mat2htmltable(tt), "</table>")
} 
if (outprn==1) {
  prn1<-FALSE; if (nx==1 & xlv[1]==0 & ny>1) prn1<-TRUE;
  tt<-c(" ",blvb); nn<-c("","Outcome","Exposure",blvb);
  for (k in (1:nclv)) {
    wdtmp0<-WD;
    if (!is.na(colvname)) {
      tt<-rbind(tt,c(clvb[k],rep(" ",nblv)))
      if (k<nclv) wdtmp0<-WD[WD[,colvname]==clv[k],];
      print(paste("Stratified by",colvname, ":", clvb[k]))
    } 
    for (i in (1:ny)) {
      if (!prn1 & ny>1) tt<-rbind(tt,c(yb[i],rep(" ",nblv)))
      for (j in (1:nx)) {
          if (!prn1) {colj<-xxlbl_[[j]];} else {colj<-yb[i];}
          nnj <- c(clvb_[k],yb[i],xb[j])
          for (z in (1:nblv)) {
            print(paste("Stratified by",bvar, ":", blvb[z]))
            if (z<nblv) {wdtmp1<-wdtmp0[wdtmp0[,bvar]==blv[z],];} else {wdtmp1<-wdtmp0;}
            wdtmp<-removeNA(i,j,wdtmp1)
            fml<-paste(yvname[i],"~",fml0[j]);
            if (z==nblv) fml<-paste(fml,"+factor(",bvar,")",sep="")
            if (k==nclv & !is.na(colvname)) fml<-paste(fml,"+factor(",colvname,")",sep="")
            if (opt=="gam") tmp.mdl<-setgam(fml,i)
            if (opt=="gee") tmp.mdl<-setgee(fml,i)
            if (opt=="glm") tmp.mdl<-setglm(fml,i)
            tmpooi<-mdl2oo(tmp.mdl,xxname_[[j]],opt)
            colj<-cbind(colj,tmpooi[[1]])
            nnj<-c(nnj,tmpooi[[2]])
          }
          tt<-rbind(tt,colj); nn<-rbind(nn,nnj)
      }
    }
  }
  if (is.na(colvname)) {nn<-nn[,-1];} else {nn[1,1]<-vlabelV[vnameV==colvname];}
  w<-c(w,"<table border=3>", mat2htmltable(tt), "</table>")
} 
if (outprn==2) {
  nn<-c(" "," ",yb);
  for (j in 1:nx) {
    tt<-c(xb[j],xxlbl_[[j]][-1]); ncum<-length(xxlbl_[[j]])-1; nnj<-xb[j]
    for (k in (1:nclv)) {
      wdtmp0<-WD; nnjk<-" "
      if (!is.na(colvname)) {
        tt<-rbind(tt,c(clvb[k],rep(" ",ncum))); nnjk<-clvb_[k]
        if (k<nclv) wdtmp0<-WD[WD[,colvname]==clv[k],];
        print(paste("Stratified by",colvname, ":", clvb[k]))
      } 
      tt<-rbind(tt,c("N(%)",varstats(wdtmp0[,xvname[j]],xxlvl_[[j]],dec)[-1]))
      for (i in (1:ny)) {
         wdtmp<-removeNA(i,j,wdtmp0)
         fml<-paste(yvname[i],"~",fml0[j]);
         if (k==nclv & !is.na(colvname)) fml<-paste(fml,"+factor(",colvname,")",sep="")
         if (opt=="gam") tmp.mdl<-setgam(fml,i)
         if (opt=="gee") tmp.mdl<-setgee(fml,i)
         if (opt=="glm") tmp.mdl<-setglm(fml,i)
         tmpooi<-mdl2oo(tmp.mdl,xxname_[[j]],opt)
         tt<-rbind(tt,c(yb[i],tmpooi[[1]][-1]))
         nnjk<-c(nnjk,tmpooi[[2]])
      }
      nn<-cbind(nn,c(nnj,nnjk))
    }
    w<-c(w,"</br><table border=3>", mat2htmltable(tt), "</table>")
  }
  if (is.na(colvname)) {nn<-nn[-2,];} else {nn[2,1]<-vlabelV[vnameV==colvname];}
}
if (outprn==3) {
  tt<-c(" ",xb[j]);  nn<-c(" ","Outcome",xb)
  for (k in (1:nclv)) {
    wdtmp0<-WD; 
    if (!is.na(colvname)) {
      tt<-rbind(tt,c(clvb[k],rep(" ",nx))); 
      if (k<nclv) wdtmp0<-WD[WD[,colvname]==clv[k],];
      print(paste("Stratified by",colvname, ":", clvb[k]))
    } 
    for (i in 1:ny) {
      coli<-yb[i]; nni<-c(clvb_[k],yb[i])
      for (j in 1:nx) {
         wdtmp<-removeNA(i,j,wdtmp0)
         fml<-paste(yvname[i],"~",fml0[j]);
         if (k==nclv & !is.na(colvname)) fml<-paste(fml,"+factor(",colvname,")",sep="")
         if (opt=="gam") tmp.mdl<-setgam(fml,i)
         if (opt=="gee") tmp.mdl<-setgee(fml,i)
         if (opt=="glm") tmp.mdl<-setglm(fml,i)
         tmpooi<-mdl2oo(tmp.mdl,xxname_[[j]],opt)
         coli<-c(coli,tmpooi[[1]])
         nni<-c(nni,tmpooi[[2]])
      }
      tt<-rbind(tt,coli); nn<-rbind(nn,nni)
    }
  }
  if (is.na(colvname)) {nn<-nn[,-1];} else {nn[1,1]<-vlabelV[vnameV==colvname];}
  w<-c(w,"</br><table border=3>", mat2htmltable(tt), "</table>")
}
sink()
prnopt<-c("β (95%CI) Pvalue / OR (95%CI) Pvalue", "β (95%CI) Pvalue", "β+se / OR (95%CI) *P<0.05 **P<0.01 ***P<0.001")

 

w<-c(w,"</br>表中数据：", prnopt[colprn])
w<-c(w,paste(c("</br>结果变量:",paste(yb,collapse="; ")),collapse=" "))
w<-c(w,paste(c("</br>暴露变量:",paste(xb,collapse="; ")),collapse=" "))
if (length(avname)==0) avb1<-"None";
w<-c(w,paste(c("</br>调整变量:",paste(avb1,collapse="; ")),collapse=" "))
if (smoothav>0) w<-c(w,". Generalized additive models were applied")
if (opt=="gee") w<-c(w, paste("</br>Generalized estimate equation were used, subject ID=", subjvname, "(", gee.TYPE,")",sep=""))
w<-c(w,"</br></br>各模型所用的样本量</br><table border=3>", mat2htmltable(nn), "</table>")
w<-c(w,wd.subset)
w<-c(w,paste("</br>此表用易侕统计软件 (www.empowerstats.com) 和R软件生成，生成日期：",Sys.Date()))
w<-c(w,"</body></html>")
fileConn<-file(paste(ofname,".htm",sep="")); writeLines(w, fileConn)